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Abstract 
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^ ' Detection of magnetic-type (_B-type) polarization in the Cosmic Microwave Background (CMB) radiation 

plays a crucial role in probing the rehc gravitational wave (RGW) background. In this paper, we propose a 

new method to deconstruct a polarization map on an incomplete sky in real space into purely electric and 

^^ ■ magnetic polarization type maps, ^(7) and ^(7), respectively. The main properties of our approach are as 

r ^ ' follows: Firstly, the fields Siff) and Biff) are constructed in real space with a minimal loss of information. 

(~| I This loss of information arises due to the removal of a narrow edge of the constructed map in order to 

I , remove various numerical errors, including those arising from finite pixel size. Secondly, this method is 

o 

;_( ■ fast and can be efficiently applied to high resolution maps due to the use of the fast spherical harmonics 

^ ' transformation. Thirdly, the constructed fields, £^(7) and B{j), are scalar fields. For this reason various 

techniques developed to deal with temperature anisotropy maps can be directly applied to analyze these 
fields. As a concrete example, we construct and analyze an unbiased estimator for the power spectrum of 
- ■ the S-mode of polarization Cf^. Basing our results on the performance of this estimator, we discuss the 

CN ■ RGW detection ability of two future ground-based CMB experiments, QUIET and POLARBEAR. 

T-H ■ 
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INTRODUCTION 



The extreme conditions in the very early Universe produce primordial perturbations of two 
generic types, namely, density perturbations (scalar perturbations) and relic gravitational waves 
(tensor perturbations) l|, y]. In the simplest scenarios, these perturbations are characterized 
by a nearly scale invariant primordial power spectra. The experimental determination of the 
parameters specifying these power spectra provides an important method to investigating the 
physics of the very early Universe. The Cosmic Microwave Background Radiation (CMB) has 
proved to be a valuable tool in this respect. Scalar and tensor perturbations leave an observable 
imprint in the temperature and polarization anisotropies of the CMB. The recent experimental 
effort, including the WMAP satellite ^, QUaD 11, BICEP ^ and so on [61-131 [d, M > lias lead 
to a robust determination of the parameters characterizing the primordial density perturbations. 
On the other hand, the detection of relic gravitational waves (RGWs) remains an outstanding 
experimental challenge, and a key task for current, upcoming and future planned CMB observations 



on the ground 
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17l |. on balloons jl8l42Cll| and in the space 211425 1. 



Both density perturbations and RGWs contribute to the various CMB anisotropy power spectra, 
namely TT (temperature), EE ("electric" -type polarization) and TE (temperature-polarization 



cross correlation) 



2614351]. In addition, RGWs produce the "magnetic" -type (BB) polarization 



that is not produced by density perturbations 



3a]. In principle, all of these information channels 



{TT, TE, EE and BB) shoul be used to infer the RGW signal in the CMB. However, if the 
contribution of RGWs to the CMB is small (r < 0.05), the BB channel will the best venue for 
detecting RGWs |35[. 

The separation of the polarization field into electric and magnetic components is a subtle issue. 
In practice, in a CMB experiment, one directly observes only the Stokes polarization components Q 
and U. Given a full sky map of the components Q and U, one can construct the so-called E-mode 
and B-mode of polarization (sometimes referred to as G and C modes, respectively) using spin- 
spherical harmonics expansion in an unambiguous manner. It is important to keep in mind that, 
by virtue of construction, E and B modes of polarization are non-local quantities, and require the 
information on Stokes parameters on the complete sky. However, in realistic cases (ground-based, 
balloon-borne and satellite experiments), the Stokes parameters are measured only on a fractional 
portion of the sky. In this situation, the simplest method for constructing electric and magnetic 
polarization fields, using the spin-spherical harmonics leads to mutual contamination, often referred 
to as the EB-mixing. This EB-mixing can become a dominant hinderance for detecting the RGW 



signal [37[. In order to overcome this difficulty, numerous methods have been developed to separate 
E and B types of polarization on an incomplete sky |38h44|| . However, these methods suffer from 



one or several of the following drawbacks - they are slow in practice, they are difficult to realize in 
pixel space, and/or they lead to partial information loss. 

In the current work, we develop a novel method to separate the electric and magnetic compo- 
nents of polarization on a partial sky. In contrast to previous works 38|-|40(] , in this paper we focus 



mainly on the construction of pure E and B types of polarization in the real space, as opposed 
to constructions in harmonics space. Our method is based on a simple redefinition of electric and 
magnetic components of polarization, so as to make them only depending on the differential of the 
Q and U fields. On a full sky, this definition is equivalent to the standard definition. However, due 
to the differential nature, our definition is directly extendable to polarization maps given on an 
incomplete sky. The main advantages of our method are as follows. Firstly, the loss of information 
is small, and is only caused by the removal of a narrow edge around the observed portion of the 
sky to reduce numerical errors. Secondly, this method is easy to realize for pixelized polarization 
maps, and is sufficiently fast so as to be practical for high resolution full sky surveys. Thirdly, 
the method leads to construction of scalar E and B type fields. For this reason, one can directly 
apply the various techniques developed for temperature anisotropy. In particular, along this route, 
we construct an unbiased estimator for the i?-mode of polarization power spectrum and gauge its 
performance in detecting RGWs. 

The outline of the paper is as follows. In Sec. m we introduce the basic notations and explain 
the general ideas behind our method. In this section we derive the basic equations, and construct 
the pure electric and magnetic fields f (7) and ^(7), given the polarization map ((5(7), ^(7)) 
on an incomplete sky. In Sec. IIIH we apply this method in practice for a pixelized polarization 
map. In this section we explicitly construct the 6(7) field, and discuss the various sources of 
contamination. We go on to discuss the effect edge removal, size of pixels and experimental beam 
size on the resulting contaminations. In Sec. IIVI we focus on the applying our method to small sky 
surveys. Combining our method with the pseudo-C^ estimator method we construct an unbiased 
estimator for the i?-mode power spectrum. We show that, in practice, our estimator performs only 
slightly worse than an estimator in an idealized situation with no loss of information. Based on this 
estimator, we analyze the ability to detect relic gravitational waves through their signature in B- 
mode of polarization in two planned ground-based CMB experiments, QUIET and POLARBEAR. 
We conclude in Sec. |V] with a brief summary of our main results. 



II. DECOMPOSITION OF THE POLARIZATION FIELD INTO ELECTRIC AND MAG- 
NETIC COMPONENTS ON AN INCOMPLETE SKY 

Let us first give a brief recount of the standard procedure to construct electric i?-inode and mag- 
netic B-mode of polarization, given a complete sky. For mathematical simplicity, it is convenient 
to introduce the complex conjugate polarization fields P± as follows 

P±il)^Q{7)±iUij), (1) 

where 7 denotes the position on the sky, and Q and U are assumed to be real fields on the sky. 
jeing ±2 spin-weighted quantities, can be expanded over appropriate spin-harmonic 



The fields P±, 

functions (see 45| for instance): 

^±(7) = X] "±2,£m ±2y£m(7), (2) 

£m 

where ±2^m(7) are the spin- weighted spherical harmonics. The multipole coefficients a±2,em are 
calculable by 

a±2,im = / d^P±{l) ±2^/^(7)- (3) 

The E and B mode multipoles are defined in terms of the coefficients a±2,e.m in the following 
manner 

Elm — ~9 [(^2,(.m + a-2,£m] , Bim = ~Tr- ['^2,em — 0.^2,lm] ■ (4) 

One can now define the electric polarization sky map E{'j), and the magnetic polarization sky map 

B{'j) as 

Eij) = Y,EemYem{7), B{^) = Y,B,mYira{l)- (5) 

Hm £m 

The power spectra of E and B modes of polarization are defined, in terms of the multipole coeffi- 
cients Elm and Bim, as 

Q = 2£ I ]_ / A^^mEim), Q = o^ ^ y .{BlmBim): (6) 

m m 

where the brackets denote the average over all realizations. 

It is important to note that, the polarization sky maps £'(7) and -6(7), are constructed out 
of underlying Q{^) and U{'y) maps in a non local manner. This is to say that, the value of the 
E 01 B field at a given point 7, in virtue of ([5]), depend on multipole coefficients E^m and -B^mj 



respectively. These coefficients, in turn, depend on integral of P±{'y) over the full sky (see ([3]) and 
dH). Therefore, one requires the knowledge of Q{^) and C/(7) (or equivalently P±{'y)) over the 
entire sky in order to construct the £'(7) and B{'j) fields. 

As was mentioned previously, in realistic scenarios, one does not have information on Q and U 
fields on the entire sky. For this reason ([SD, (HD and (0) cannot be applied directly to construct 
E and B types of polarization maps on an incomplete sky. In order to avoid this problem, in the 
present paper we adopt a different but related definition for electric and magnetic polarization 
maps 
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-[gig2P+(7) + SiS2P-(7)], 
-i[gig2P+(7)-SiS2P-(7)], 



(7) 
(8) 



where g^ and ds (s = 1, 2) are the spin lowering and raising operators, respectively 



/ ■ /ix-s I ^ id 



sm 



[sin" 6] 



(9) 
(10) 
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and ([8]) for £ and B have been previously discussed in the literature 4a, 



481 ]. These have often been denoted as E and B, respectively (see for example 



48l|). The fields £ and B are equivalent to the fields E and B introduced in [30], where 



these fields were introduced as two independent invariants constructed out of the second covariant 



derivatives polarization tensor (see Eq. (36) in [30|). In the present paper, to maintain a clear 



distinction from {E, B) in ([5]), we shall use the (£", B) notation. 

The constructed electric and magnetic fields are scalar fields on the sphere. Therefore, assuming 
£ and B given on a full sky, one can determine the spherical harmonics decomposition coefficients 



£,^ ^ / d^ £(j)Yl^{i), Berr, = / dj ^(7)5^4(7) 



(11) 



These relations can be inverted to give 



^(7) = X] ^^™^^™(7)' ^(7) = ^ BemYemi'l)- 



(12) 



Im 



Im 



The multipole coefficients £im and Bim are related to Eim and Bim defined in (jlj) by a ^-dependent 
numerical factor Ng = y/(^ + 2)!/(^ - 2)! [ad, M, m 



£em = NeEim, Bem = NeBi„ 



(13) 



One can also define the power spectra of 8 and B modes of polarization 



m m 

These are related with the power spectra Cf and Cf through relations 

Cl' = NJCf'^, c!^ = Nfcr- (15) 

Thus, in comparison with Cf^ and Cf^, the power spectra Cf and Cf^ are "bluer", due to the 
factor Nf. Note that, the relations (J13p and (jlSp assume the polarization fields to be given on a 
complete sky. 

It is important to point out that, the quantities 8 and B defined in ([7]) and ^ only depend 
on the differential of the Q and U fields by construction. Therefore, these definitions can be, in 
principle, applied in the case of Q and U given on an incomplete portion of the sky, to construct the 
8 and B fields on this portion. We now proceed to discuss the relevant steps for this construction 
on an incomplete sky. 

In order to describe the partial sky observations, we firstly introduce the mask window function 
W{'^). This mask function is non-zero only in the observational region of the sky. In addition, we 
shall assume that VF(7) is a real function. In particular, the special case with VF(7) = 1 in the 
observational region corresponds to the widely discussed top-hat window function. In the present 
paper, we denote this special case of a top-hat window function as w{'y). With the introduction of 
the window function W{'^), the analysis of the polarization field P±{'^) defined on the partial region 
of the sky becomes equivalent to studying the masked field P±{'))W{'^) defined on the complete 
sky. 

In the general case of an arbitrary mask, one can define two full sky maps 8{'y) and Bify) 
constructed out of observational data 

m ^ -\ ^i^2{p^(ci)wm + gig2(P-(7)w^(7))] , (16) 

^(7) = -\. [8i82(P+(7)VF(7)) - ^x^2{P-{n)W(ci))\ . (17) 

Due to the presence of the window function W{f^\ the two maps 8{pj') and B{^^ do not correspond 
to pure electric and magnetic types of polarization. The main task of this work is to construct 
pure £^(7) and ^6(7) fields out of f (7) and ^6(7). 

Moving on, we define the multipole decomposition coefficients Eim and B^^ as 

^£m = ^ I ^7 m ^4(7), ^^rn = ^J ^7 ^(7) YL^l)- (18) 



With this definition, the E and B fields can be expanded in terms of the multipole coefficients E'^m 
and Bini in the foUowing manner 

m) = Y.N,E,^Y,^{^), B{^) = Y,NiBimYUl)- (19) 

im £m 

The multipole decomposition coefficients Egm and Bim can be calculated in an alternative manner. 
One can begin by defining the complex polarization fields P±{j) = ^±(7)^^(7) and construct the 
multiple coefficients using ([3]) and @ (with tilde placed on all the relevant quantities). It can be 
verified that the two definitions are equivalent. 

Before proceeding, let us point out some simplifying relations. Firstly, from the definitions of 
Bs and Qs, it follows that §2 = (S2)* and §1 = (9i)*. In light of definitions ([7]), ([8]), ([H]) and (fTT]) . 
it follows that the two sets of fields {£, B) and (£', B) are real. Furthermore, from the definitions 
one has 

8182-P+ = -£- iB, diB2{P+W) = -£- iB. (20) 

Thus, in order to determine the relation between the two sets of fields (£', B) and (£", B)^ it suffices 
to study the relation between §182^+ and 8i82(-P+VF). 

In order to derive the relation between the two sets, it is convenient to expand the quantity 
[8i82(i^+W^)]VF, using the definition of 8^ operator ^ in the following way 

[^i^2{P+W)]W = {^i^2P+)W'^ + {^iW){^2P+)W + {d2W){^iP+)W 
+ P+W{^i^2W) + cot ^[1^(82^+) + P+{^2W)]W 
+ {2 + 2 coi^ 9) P+W'^ + 2 coieWdi{P+W). (21) 

Using the following set of relations that follow from ^ 

(82^+)^ = ^2{P+w) - P+82VF - 2 cot ep+w, 
{^iP+)w = ^i{p+w) - p+diW - cot ep+w, 

^^(82^+) + P+{^2W) = 82(P+W") -2coieP+W, 

along with expressions in (|20p we arrive at the expression 

[£ + iB]W'^ = [£ + iB\W + i{^iW)[d2{P+W) - P+^2W -2 cot 6 P+W] 

+ {d2W)[^i{P+W) - P+^iW - cot OP+W] + P+W{^i^2W) 
+ coie^2{P+W) + 2W'^P+ + 2Wcoie^i{P+W)\ . 

(22) 



This expression can be rewritten in a compact form 

[£ + iB]W'^ = [£ + iB]W + ct, (23) 

where ct denotes the correction term. This correction term is complex in general. The real and 
imaginary parts of the correction term are given as 

Re[ct] = Q[3cot OWW^ + W{W^^ - Wyy) - 2{W^ - W^)] 
+ U[2 cot eWWy + 2WW^y - 4W^Wy] 
+ 2W^ \{QW), + {UW)y\ + 2Wy[{UW)., - {QW)y] , (24) 



and 



Im[ct] = U[3 cot eWWa; + W{W^^ - Wyy) - 2{W^ - W^)] 

- Q[2 cot eWWy + 2WW^y - 4Wa;Wy] 

- 2Wy[iQW), + {UW)y] + 2W,[{UW),-{QW)y]. (25) 
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In the above expressions we have introduced the shorthand notations F^ = ^, Fy = -^ 
Fxx = ^, Fyy = ^. Jgg 2 and F^y = ^-JqLqq for an arbitrary function ^(7). In Appendix [Al we 



sin 

discuss the question of numerically calculating the various terms in the above expression in pixel 
space. 

Finally, one can construct the pure electric and magnetic fields £ and B on the observed portion 
of the sky (i.e. region of the sky for which 1^(7) 7^ 0) using expression 

[£ + iB] = [£ + iB]W-^ + ctT^-2. (26) 

The construction of the pure electric and magnetic scalar fields £ and B is the main result of this 
paper. It is worth pointing out that the construction of these fields is independent of the choice of 
the mask function W{^), as long as the mask in non-zero in the observed portion of the sky. This 
method for recovering the scalar fields £ and B is lossless in the real space in following sense. If one 
was given the polarization fields Q and U on the entire sphere and constructed the corresponding 
B field using ([8]) (or the £ field using ([7]) ) and compared the resulting scalar field in the observed 
region with result of above procedure (j26p . one would find the two fields equal. However, due to 
the ill-behaved nature of W~^ and W~'^ at the edge of observed region, it is difficult to realize the 
above construction in practice. In order to circumvent this problem, as will be discussed in the 
following section, one has to remove the edge of the constructed polarization maps. 



In conclusion of this section it is instructive to clarify the issues associated with possible leakage 
from the so-called ambiguous modes. It is known that, on a manifold with a boundary, the decom- 
position of the polarization field, in addition to pure E and B components, contains ambiguous 
modes that satisfy both E-uiode and B-mode conditions simultaneously (see 



381, 



39[ for details). 



In particular, when constructing the power spectrum estimators for B mode one has to ensure that 
there is no leakage from the ambiguous modes. In the current work, the Bj^eci'j) field does not 
contain contribution from either E modes of polarization or ambiguous modes, by the virtue of 
construction (analogous to XB in [41]). For this reason, the power spectral estimators constructed 
from this field will be free from contaminations from both i?-modes and ambiguous modes. 

III. £/B SEPARATION IN PIXEL SPACE 

In this section, we shall discuss the issues related to separation of electric and magnetic polar- 
izations <?(7) and 6(7) in the pixel space using the results of the previous section, in particular 
expression (p6|) . We shall discuss this procedure using a toy model. For this toy model, we assume 
that an ex per iment will only observe the Stokes parameters Q and U in the northern hemisphere. 
Following [371, |40| , we adopt the following axially-symmetric form for the mask window function 
1^(7) 



W{j) 



( 

1 




0< 00-^1, 


1 

2 


-icos(^vr) 


00 - ^1 < ^ < 00, 







0>0o. 



(27) 



In the above expression, 0o corresponds to edge of the observational area, and 0i > is the 
smoothing scale. The limiting case, 0\ = 0, corresponds to the top-hat window function (ii'(7) = 1 
for 6 < Oq and w{^) = for 6 > 0q). However, the top-hat function is discontinuous at 6 = 0q, which 
makes quantities Wxx, Wxy and Wyy ill-defined at the edge. In order to avoid these difficulties, 
it is convenient to use a window function with 61 ^ 0. Throughout the present section we use 
values 9q = 90° and 9i = 30°. It is important to point out that the formalism outlined in Sec. |II] 
is applicable for arbitrary mask window functions, not necessarily axially symmetric. The simple 
symmetric form (|27p for the mask was chosen for simplicity and clarity of presentation. In realistic 
scenarios one will have to use a more complex mask, that will take into account the non-symmetric 
form of the observed region and various point source contaminations. 

For an axially symmetric window function 1^(7) (i.e. when 1^(7) is independent of </)), such as 



the one considered in the present section, the correction term ct in ()24p and (j25p simpUfy to 

Re[ct] = Q[3coteWW^ + WW^^-2W^] + 2W^[{QW).^ + {UW)y], (28) 

Im[ct] = U[3 cot eWWcc + WW^^-2W^] + 2W^[{UW)x-{QW)y]. (29) 

In order to demonstrate the £/B separation we shall work with simulated polarization maps. 
We use the synfast subroutine included in the HEALPix package to generate a full sky map of (3(7) 
and U{'y) fields. In order to generate this map we use the best-fit WMAP5 values for cosmological 



parameters 



49| 



Qbh'^ = 0.02267, Qch'^ = 0.1131, Qa = 0.726, 



Treion = 0.084, A^ = 2.446 X 10"^ n, = 0.96, (30) 

and assume no contribution from gravitational waves and cosmic lensing to the B-mode of po- 
larization, i.e. C^^ = 0. We adopt the pixelization with Nside = 512. We set the full width at 
half maximum (FWHM) for the Gaussian beam to Op = 30'. In what follows we shall be solely 
interested in the determination of B polarization field, and studying the possible contaminations 
to it. In this context, since Cf^ = 0, a non-zero field B would be wholly attributed to the residual 
£^i3-mixing contamination. 

Before proceeding to construct the pure magnetic field B, we shall construct and analyze the 
auxiliary B field. From the simulated Q{^) and U{^) maps we construct the ^6(7) in the following 
manner. Firstly, we construct the multipole coefficients Bim- This is done by building multipole 
coefficients d±2/m = J dj -P± (7)^^(7) ±2^£m(T') using the simulated ((5(7), ^^(7)) maps along with 
the window function W{j) given in (p7|) . and then calculating Bi^ — ~Yi [^2/m — a-2,^m]- These 
steps were performed numerically using the anafast routine from the HEALPix package. Following 
this, we use (|19p to construct B from the multipole coefficients -B^m- The resulting B map is 
illustrated in Fig. [TJ 

One can see that although the map was generated with Cf = 0, B y^ in the region 9q — 6i < 
9 < 6q (i.e. 60° < 9 < 90°). This can be viewed as a of leakage of iS-type of polarization into 
B field, due to the presence of the window function W. In order to quantify this leakage in the 
harmonic (multipole) space, we construct the pseudo power spectrum as 



^' =2£ + 



^ Y, B^mBlm^ where B,^ ^ f mY*U^)dj. (31) 

m 

The resulting pseudo power spectrum is plotted a the black line in Fig. [3j 
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FIG. 1: The B map constructed from an input model with no magnetic polarization (in /iK). 
We can now reconstruct pure magnetic type field B using (|26|) in the following manner 

^rcc(7) = m)W-H^) + Im(ct(7))Ty-2(7). (32) 

We use the subscript rec to indicate that this field was reconstructed field from B and the {Q, U) 
fields. The results of reconstruction are presented in Fig. [2j Since the input cosmological model 
assumes no contribution from gravitational waves, one can expect B^ecij) = 0. A visual comparison 
of Fig. [Hand Fig. [2] shows that we have been able to remove much of the leakage that was present 
in B. The remaining residual contamination in B^cc is shown (with a magnified scale) in the right 
panel of Fig. [2j This residuals are a small fraction (~ 2%) of the total leakage in Fig. [TJ In order to 
quantify these residuals in multipole space, we construct the pseudo spectral estimators replacing B 
with Brec in ([3T]) . The resulting pseudo power spectrum is plotted as the red curve in Fig. [3l It can 
be seen that the resulting leakage power for Bj^ec is significantly lower that the corresponding power 
for B. In particular, in the practically interesting range of multipoles £ G (50, 200), the spectrum 
for the reconstructed B^ec field is roughly four order of magnitude lower than the spectrum for B. 
The remaining residuals in ,Brec are attributed to numerical errors. 

We believe that the small remaining residuals in B^cc are a result of two types of numerical 



39|, it was 



errors that cannot be avoided in practice. The first reason for errors is pixelization. In 
argued that pixelization can lead to the mixing of electric and magnetic modes. This point can be 
intuitively understood in the following way. Imagine a survey that observes polarization on a small 
square region of the sky. Pixelization introduces a Nyquist wavenumber kj^, such that all modes 
with wavenumbers greater than k^ are aliased to modes with wavenumbers less than the Nyquist 
value. This aliasing completely shuffles the direction of wavenumbers, thus essentially leading to 
a complete mixing of electric and magnetic modes. Although the complete avoidance of the errors 
due to pixelization is impossible, these numerical errors can be reduced by using a larger value for 
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FIG. 2: The pure magnetic field ,Si.cc(7)- The edge of the map with W{'j) < 0.03 has been removed. The 
left panel is scaled similarly to Fig. [TJ while the right panel has the scaling magnified in order to show the 
residual leakage. Both the panels use the units /iK. 



-^side (see Sec. IIII Bl for details). In the present evaluation, with A'side = 512, finite pixelization 
seems to be the main reason for residual power spectrum for i > 150. The main contribution to 
this residual power spectrum comes from the residual Brecij) at the pole in the real space (see the 
right panel in Figure [2]). 

The second reason for numerical errors is the steep growth of power spectrum Cf with in- 
creasing i. Due to this, even a small relative numerical error at higher multipoles seeps through to 
lower multipoles. In other words, these errors occur due to fact that the various sources of noise 
and E mode signal are not band limited. We believe that these types of errors mainly account for 
the residual power spectrum of Brcc at multipoles I < 150. These errors are mainly caused by the 
residual Bj-ecil) around the observed edge 9 = 9q m the real space (see the right panel in Figure 

ED. 



A. Dependence of residual leakage on edge removal 



Figure [2] shows that much of the residual leakage occurs around the edge of the observation 
region. It is therefore instructive to study the edge effects in more detail. The expression (j32p for 
;Brcc depends on Im(ct) correction term given by pop . This correction term contains two terms, 
(3 cot eWWx + WW^^ - 2Wl)lW'^ and IW^jW'^, which have the following asymptotic at the edge 
of the map as ^ — t- ^q 

'izotQWW^^WW^^-lWl 6 2W, 160? 

Thus, the two functions are divergent for 6 —^ 9q. This implies that the signal-to-noise will tend 
to zero for data as the boundary of the observed region is approached. Because of this divergence, 

12 




J I I I I ' ' I 



_l I I I I I I— L. 



1000 



Multipole: I 



FIG. 3: The pseudo power spectrum Df ^ (black line, i.e. line 1) of the B field and pseudo power spectrum 
Di^ (red line, i.e. line 2) of the Bicc field. 

in numerical calculations, one must remove the edge of the map in order to avoid numerical errors 
associated with these divergences, thereby introducing a small loss in information. 

In order to investigate the dependence of residual leakage on the edge removal, in Fig. H] we 
plot the pseudo power spectrum D^ of residual Brec constructed for the same simulated data but 
for two different edge removals. The first case (red line) corresponds to portion of the sky with 
W < 0.03 removed (this corresponds to removal of data with Oq — 9 < 0.06). The second case 
(green line) corresponds to a portion of sky with W < 0.1 removed (corresponding to removal of 
data with ^o ~ ^ < 0.1). The second case corresponds to a larger portion of the sky removed, and 
thus to a larger loss in information. This loss of information leads to a smaller value for the power 
spectrum at lower multipoles I < 150. This fact can be clearly seen from right panel of Fig. [2l As 
one removes more of the data from the equator the residual spectrum of B^ec becomes smaller. On 
the other hand, in the region of higher multipoles, where the dominant contribution comes from 
finite pixelization errors, the two power spectra are comparable. 
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Multipole: I 
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FIG. 4: The red line (i.e. line 2) shows the pseudo power spectrum Dg calculated for the edge removal at 
W — 0.03, while the green line (i.e. line 3) shows the De for edge removal at M^ = 0.1. For comparison, the 
black line (i.e. line 1) shows the pseudo power spectrum Df^. Note that, the black line (i.e. line 1) and red 
line (i.e. line 2) are identical to those in Fig. [31 



B. Dependence of residual leakage on pixelization number TVgi 



idc 



As was pointed out earlier, one of the reasons for residual leakage of power into B^cc is finite 
pixelization of the sky map. In order to demonstrate the effect of pixelization on the residual 
power, in Fig. [5] we show the pseudo power spectrum D£ calculated for two different pixelization 
numbers A'side = 512 (red line) and A^side = 1024 (green line). As one might expect, the increase 
in the pixelization number reduces the leakage power spectrum. This reduction is most dramatic 
at higher multipoles i > 150, where it is two orders of magnitude for this example. At lower 
multipoles £ < 150 the reduction is not as dramatic, and is roughly by a factor 3. These results 
are consistent with our previous statements about the cause of numerical errors. Indeed, at higher 
multipoles the main cause of errors seems to be finite pixelization, whereas at lower multipoles the 
errors are generated by a combination of factors. 
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FIG. 5: The red line (i.e. line 2) shows the pseudo power spectrum D^ calculated for A'sidc = 512, while the 
green line (i.e. line 3) shows D( calculated for A^sidc = 1024. For comparison, the black line (line 1) shows 
the pseudo power spectrum Df^. Note that, the black line (i.e. line 1) and red line (i.e. line 2) are identical 
to those in Fig. [3] 

C. Dependence of residual leakage on Op 

The full width at half maximum parameter Op has an important effect on the residual leakage 
of power into B^cc- In order to understand the reason for this, one has to remember that one of the 
two reasons for residual leakage is the steep growth of power spectrum Cf^ with increasing I. The 
parameter Op regulates the exponential damping of this power spectrum at multipoles £ c^ Op , 
and therefore limits the propagation of the power in Cf into D£. 

The various contributions to the spectrum Cf^ are illustrated in Fig. [UJ The main contribution 
to power spectrum comes from density perturbations (dashed blue line). For comparison, on this 
figure, we show the contribution to Cf^ (solid blue line) and Cf^ (solid red line) from gravitational 
waves (characterized by tensor-to-scalar ratio r = 0.1), as well as contribution to C^^ from lensing 
(red dashed line). The spectrum Cf^ from density perturbations at high multipoles acts as the 
main source for the residual leakage into Di. 
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FIG. 6: The polarization power spectra generated by density perturbations (d.p.), gravitational waves (g.w.) 
with r = 0.1, and cosmic lensing (lens). 

The contribution to the various spectra at high multipoles are effectively damped by a choice of 
an appropriate Op- This parameter leads to the damping of power spectrum Cf^ proportional to 



exp 



Un2 



In Fig. [7] we show the residual leakage for two choices of the FWHM parameter 



Op = 30' (red line) and Op = 10' (green line). For comparison, in this figure we also show the 
pseudo spectrum Cf calculated Op = 30' (black line) and Op = 10' (blue line). As one might 
expect, the residual power spectrum reduces significantly with an increase in Op. For this reason, 
for the purposes of extracting the magnetic pattern of polarization in experiments with small Op 
(for example POLARBEAR experiment discussed in Sec. lIVEj) . it becomes necessary to artificially 
increase Op in order to reduce residual leakage in Srec- In Appendix[Dl we suggest a 'map smoothing' 
technique to achieve this goal. 
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FIG. 7: The red line (i.e. line 2) shows the pseudo power spectrum Di, calculated for a map with 9f = 30', 
while the green line (i.e. line 4) shows D( calculated with Op = 10'. The black line (i.e. line 1) shows Df^ 
calculated for 6p = 30', and the blue line (i.e. line 3) shows Df^ for 9p = 10'. Note that, the black line 
(i.e. line 1) and red line (i.e. line 2) are identical to those in Fig. [31 

IV. £/B SEPARATION AND POWER SPECTRUM ESTIMATION FOR SMALL SKY 
SURVEYS 

In Sec. HI] and Sec. lIIII we developed a method to construct pure electric S^cc and magnetic fields 
Brec out of the original Stokes parameter fields Q and [/ on a fractional portion of the sky. Ignoring 
the small numerical errors, it was shown that the resulting fields did not exhibit mixing. A crucial 
point about the constructed fields is that they are scalar fields. For this reason, one can use all of 
the robust techniques developed for studying CMB temperature anisotropy to the fields S^ec and 
Brec- Based on appropriation of these techniques, in this section, we shall focus on an important 
practical application, namely constructing the estimator for the power spectrum of the i?-mode of 
polarization Cf . For this reason, as in the previous sections, we shall restrict our analysis to just 
the magnetic field B^ec- 

The question of constructing an estimator for the power spectrum Cf from the field Srec 
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is analogous to the problem of construction an estimator for the temperature anisotropy power 
spectrum Cj'^ given a temperature map on a partial sky . Fortunately, there are a large number of 
methods that have been developed for this purpose [5CH54i| . Amongst these, a popular method is 



the so-called 'pseudo-C^' estimator method [53|]. This method can be easily realized in pixel space 
using fast spherical harmonics transformation, and has been applied to various CMB observations 
including WMAP data [5.5|]. However, it is well known that pseudo-C^ estimators are sub-optimal, 
particularly for low multipoles. For this reason, many authors have developed alternative estimators 
that are optimal, in particular the maximum likelihood estimators in pixel realization 



5C 



52]. The 



fundamental problem with the maximum likelihood estimators is that these methods are very 
slow, especially for larger multipoles. For large sky surveys, such as the Planck satellite, the use of 
hybrid estimators, which combine the two methods, has been suggested 5j]. The hybrid estimator 
combines the best of two worlds, it is nearly optimal and can be realized of a laptop computer even 
for large sky surveys such as Planck. 

In the present section we shall focus on small sky polarization surveys, corresponding to various 
ground-based CMB experiments J, la, Il2l4l 



14l ]. Since these surveys will be primarily sensitive to 
relatively large multipoles (. > 20, we shall restrict our analysis to pseudo-C^ type estimators, which 
are nearly optimal for large multipoles. Before proceeding we would like to point out that a hybrid 
type of estimator could be potentially used to construct an estimator for Cf^ from 6rec in the 
case of large sky surveys such as Planck. We leave this exercise for future. 

Below we shall work with a small fraction of the sky characterized by a window function ()27p with 
6*0 = 20° and 9i = 10°, corresponding to a 3% sky survey. In an ideal case, neglecting numerical 
errors, the reconstructed field -Brec(7) would be related to the underlying full sky field B{'^) through 
-SreclT) = B{l)w{'^), where '^(7) is the corresponding top-hat window function. However, as was 
pointed out in Sec. IIII Al one needs to remove a narrow edge from the observational area in order 
to avoid excessive numerical errors. For this reason, in practice, we remove the region 9q — 6< 0.03 
(corresponding to 9 > 18°) from the analysis. Below we shall use the notation w'^^) to denote the 
top-hat window function for data with edge removal. 
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A. Pseudo estimators 

The first step in constructing the pseudo estimator is the definition of spherical harmonics 
coefficients a^^ of the scalar field Srec(7) as follows 

d7^rcc(7)>V(7)5^4(7), (33) 

where W(7) is the weight function. In principle, one can choose an arbitrary form for the weight 
function. In particular, the choice W(7) = 1 corresponds to the widely discussed pseudo-C^ 
estimator introduced in [53]. This choice will be the main focus of our attention in the present 
work. An alternative choice Wij) = 1^(7) (where 1^(7) is the mask window function in Eq. (j27p ) 
corresponds to the analysis in 141 1. where it was shown that the resulting a^m lead to the pure 
B-mode estimators defined in [40]. The comparison of this choice for the weight function with our 
main choice W(7) = 1 is discussed in Appendix [Cl Th e optimal choice of the weight function in 



various cases has been discussed in 



37 



41 



56( 1 . In [4]| the authors suggest a general method to 
build the weight function VV(7) for different multipole i in order to optimize the estimator. At this 
point it is important to emphasize that although B^cci'^) preserves the available information in real 
space, a non-optimal power spectrum estimation will lead to loss of some of this information. This 
makes the study of the optimal choice of weight function particularly important. However, in the 
current paper we concentrate mainly on the simplistic case W(7) = 1, leaving the important but 
complicated question of optimal choice of weight function for future work. 

For the choice W(7) = 1, the spherical harmonics coefficients aim in dMD take the simplified 
form 



1^76^00(7)^4(7)- (34) 



These are related to the coefficients Bim (which were defined in (jlip in terms of the underlying full 



sky map -6(7)) through the coupling matrix Kirni'm' (see for instant 



m 



aim = 2_^ BiiBi'm'Kimi'm' = 2^ B (^i N (_i B (,i m' K irn& m> ■, (35) 

I'm' I'm' 

where Bi is a window function describing the combined smoothing effects of the beam and the 
finite pixel size. The coupling matrix K can be expressed in terms of the function w'{j) as 

Ke,rme2m2 = dl w' {'^)Ye,^.,nAl)Yi2'm2(l) ■ (36) 



19 



The pseudo estimator D^ is defined analogous to (|3T|) in terms of the multipole coefficients ([3 
as 



m 



— Y 2^ aimolm- (37) 



Using relations Q, (jlSp and (j35p . one obtains that the expectation value of this estimator D^ is 



related to the true power spectrum Cf^ by the following convolution 



{Di) = Y^ Mu'Blcf,^ = Y^ Mu'NIbIC§^. (38) 

The coupling matrix M in the above expression can be expressed in terms of 3j symbols as 

M,.,,.(2.. + l)E^^»J''''M^. (39) 

where tf^ is the power spectrum of the window function w'{'^) defined in an analogous manner to 

It can be shown that the covariance matrix for the pseudo estimator D^ has the form 
i^D.^D,) = .,,^'^,. E E E BlNlCl-BlNlCf^- x 

mm' £imi £21712 

X Kgmeimi K^ijn'eimi -^^m^2m2 ^i'm'e2m2 ■ 

(40) 

As it stands, this formula is not useful due to the high cost of computation. However, for high 
multipoles, this formula simplifies to [54l | 

{ADiADe) « 25|iV|cf ^5|Ar|c|^M«//(2/ + 1). (41) 

In order to implement and verify the above analytical results we have conducted numerical 
calculations using simulated data. In the first instance, we generate 1000 random full sky (Q, 
U) maps with no contribution from gravitational waves (i.e. r = 0) and no lensing. For each 
realization, we reconstruct the magnetic field B^ecij) and evaluate the pseudo estimator Dg. The 
average over 1000 realizations Di is plotted (green line) in Fig. [8l Note that, here and below we 
use the over-line to denote averaging over simulated realizations, as opposed to the angle brackets 



which denote ensemble averaging. The average for the uncleaned spectrum Df^ (defined in pip ) 
is plotted (black line) for comparison on the same figure. For next calculations, we simulate 1000 
random full sky maps with contribution from gravitational waves characterized by r = 0.1 and 
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contribution to i?-mode of polarization from cosmic lensing. The average value of the estimator 
Di is plotted (red line) in Fig. [5J The comparison of curves in Fig. [5] shows that the residual noise 
contribution to the pseudo estimator due to numerical errors (green line) is negligible in comparison 
with the contribution to the estimator from the signal (red line). One can therefore conclude that 
the resulting pseudo estimator Dg^ is effectively free from EB-rmy.m.g. 

In Fig. m in order to verify (|38|) . we plot the left-hand side (solid blue line) of this equation 
for Cf a model with r = 0.1 and contribution from cosmic lensing 57|. For comparison, in 
this figure, we plot the average -D^ over 1000 realizations for the same model (solid red line). As 
expected the two lines are close to each other, being practically indistinguishable for multipoles 
H. > 20. For comparison, in Fig. [9l we also plot the individual contributions from gravitational 
waves (solid magenta line) and cosmic lensing (solid green line). Finally, in Fig. [9l we plot the 
square root of the average over 1000 realizations of the diagonal terms in the covariance matrix 

,1/2 



yAD(ADij (dashed red line). In order to check the analytical approximation (j4ip . we also 
plot the diagonal term (AD^AD^)^'^ evaluated using the right side of expression ()4ip (dashed blue 
line). As expected the two curves for the covariance matrix practically coincide for large multipoles 
i > 80, which corresponds to the region of applicability of the approximation (|4ip . 

B. Unbiased estimators for Cf^ 

Having constructed the pseudo estimator D^, we are one step away from constructing an un- 
biased estimator for the B-mode power spectrum Cf . In this subsection, we shall discuss this 
construction. Let us for the moment assume that the coupling matrix Mi£i in l\'6S\) is invertible. In 
this case, using relation (j38p . one can immediately verify that the estimator defined as 

is an unbiased estimator of the power spectrum Cf^. In practice, this simple estimator can be used 
in large sky surveys such as Planck, where Mw is indeed invertible. However, in the case of small 
sky surveys, one cannot construct this estimator since the matrix Mu' becomes singular. In this 
case it is possible to bin the pseudo estimator data into multipole bins, and construct an unbiased 
estimator for the binned power spectrum. Following the analysis for temperature anisotropy 58|, 
we build the so-called full-sky CMB bandpowers P^ as 

Pr = ^M,-ij;p,,,Z?,, (43) 
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FIG. 8: The averaged pseudo estimator De from 1000 simulations. The red line (i.e. line 2) shows the result 
for an input model with r = 0.1 and contribution from cosmic lensing. The green line (i.e. line 3) shows the 
result for an input model with no magnetic polarization (i.e. Cf^ = 0). For comparison, the black line (i.e. 
line 1) shows the averaged estimator Df^ for an input model with no magnetic polarization. 



where the subscript b denotes the multipole bands, pm is a binning operator in £-space defined as 

ow — 1( 

otherwise 



Pbe = < 






low low 





(44) 



The non-singular binned coupling matrix M^h' participating in (143p is constructed from the coupling 
matrix M«/ 



Mw = ^PM^Mu'B, 



ele.'b'- 



(45) 



The function Bf, takes into account the effects arising due to finite beam size and finite pixelization. 
In the above expression, qib is the reciprocal operator of pu 



qeb 



2^^l if /C^) <f< /(^+1) 

otherwise 



(46) 
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FIG. 9: All of the lines are a result of averaging over 1000 simulated samples. The red line (i.e. line 2) 
shows the averaged pseudo estimator De for an input model with r — 0.1 and contribution from cosmic 
lensing. The red dashed line (i.e. line 3) shows the square root of the average of the diagonal terms in the 

,1/2 



covariance matrix (AD^AZ?^) . The blue solid (i.e. line 1) line shows the analytical result of {Di), for a 
model with r = 0.1 and contribution from cosmic lensing. The contributions two contributions are shown 
separately with a magenta solid line (i.e. line 5) for the gravitational waves and a green solid line (i.e. line 
6) for the cosmic lensing. The blue dashed line (i.e. line 4) shows the analytical approximation (|4ip for 

{ADeADey/^ 



It is straightforward to verify that P^^ is an unbiased estimator of the i?-niode of polarization 



power spectrum £{£ + 1)C^ /lir, i.e. 



(pD 



+ 1) ^BB 



27r 



-Ci 



The covariance matrix of the bandpowers is related to the covariance matrix (AD^ADg/) in (^01 



by 



59| 



(APf^AP^^) = M^,lp,^e{ADeADe){pb,e'f{M-^f- 



(47) 



In Fig. [To] we plot the value of the bandpower Pf ^ (red dots) averaged over 1000 realizations. 
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The realizations were generated for a model including contribution from gravitational waves with 
r = 0.1 and cosmic lensing. The multipole bins were chosen with A£ = 10 for each bin. The 
analysis shows that (up to discrepancies that can be attributed to finite number of realizations) 
the average of the power spectrum estimators coincide with the theoretical (input) spectrum. The 

/ X 1/2 

error bars ( AP^^AP^-^ I (red error bars) were calculated using (I47p . with ensemble average 



b '-^'^ b' J 
replaced by an average over realizations. As can be expected, the error bars are large for the 

first three data points, due to the small sky coverage. In addition to evaluating error bars, it was 

verified that the correlation between various multipole bins is QUIET weak (all of the correlation 

coefficients are smaller than 0.3). Note that, the correlation matrices and corresponding error bars 

calculated here do not include contribution from instrumental and astrophysical foreground noises. 

The left panel in Fig. [10] shows the power spectrum estimation with both gravitational wave and 

cosmic lensing contributions included. The right panel shows the power spectrum estimation for 

the gravitational wave contribution alone. In the right panel, cosmic lensing serves as an effective 

noise for the detection of gravitational waves (see next subsection for details). It can be seen 

that the gravitational wave signal is larger than the corresponding error bars only in a range of 

multipoles 50 < ^ < 150, peaked at ^ ~ 100, consistent with results in [60]. 

In order to quantify the detectability of gravitational wave signal, it is convenient to introduce 

the total signal-to-noise ratio as follows 



^/^ = ./E<Pr(g^))(Cov-^)66'(Pr(gw)), (48) 

V bb' 

where Covf,;,' = (AP^^AP^^) is the covariance matrix of the bandpower estimator (j47p . For the 
example considered above with r = 0.1, we find S/N = 8.26. 

It is important to emphasize that our pseudo-Cf estimator is QUIET different from the pseudo- 
CV polarization estimators suggested in [37|, 



m 



61[, or an equivalent estimator suggested in (gj]. In 
6ll ]. the unbiased estimators are constructed directly from the pseudo estimators Cf and 
Cf , both of which are a mixture of electric and magnetic types of polarization. The resultant 
mixing increases the magnitude of the covariance matrix for the unbiased estimator, and becomes 
one of the main contaminations for the detection of gravitational waves. In [371], the authors found 
that, for small sky surveys covering one or two percent of the sky, the mixing contamination to 
the covariance matrix of B-uiode power spectrum estimator typically limits the tensor-to-scalar 
ratio that can be probed to r > 0.05. On the other hand, the pseudo-C^ method suggested in the 
present work explicitly separates the electric and magnetic types of polarization up to very small 
numerical errors. For this reason, the effects of mixing of electric and magnetic modes, which are 
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completely removed (reduced to negligible levels) in our case, do not put a limit on the ability to 
detect gravitational waves. This is the main advantage of our method, and is the main motivation 
for this paper. 

At the end of the subsection, we would like to point out that, if one proceeds to construct 
an unbiased estimator for Cf^ using the Brec{'y)W {'y) field (instead of ;Brec(7) field used above), 
the resulting unbiased estimator will be equivalent to the "pure B-mode" estimator defined in 
[40|. This has been discussed in Appendix [Cj On the other hand, if one constructs the unbiased 

estimator for C^^ using 0(7) (instead of Brccil)) adopting a top-hat window function (instead of 

i I , , 

W{'y)), one would return to the i?-mode estimator defined in f44]. As was pointed out in [4^, the 

resulting estimator suffers from large EB-mixing at the edge of the observed field. 

C. Information loss due to edge removal 

As was emphasized in Sec. IIIIl in practical calculations, the edge of the partial sky map has 
to be removed in order to reduce numerical errors. The edge removal leads to the partial loss of 
information. In this subsection we study the impact of this information loss on the performance of 
-B-mode of polarization power spectrum estimator P^^. 

In order to study the performance of the estimator P^^ in an ideal case with no edge removal 
we perform the following steps: 

1. We generate 1000 full sky (Q, U) maps, for a cosmological model with r = 0.1 and contri- 
bution from cosmic lensing. For each of these maps we calculate the multipole coefficients 
Be^ using ([ID, ([SD and ©. 

2. With the multipole coefficients B^m we construct the full sky map 6(7) using (I12p . 

3. We construct the top-hat mask window function ^(7) equal to unity for 9 < 9q = 20° and 
zero otherwise. We now construct the masked magnetic field Brecil) = ^^(7)^(7). The 
masked field ;Brec(7) constructed in this manner corresponds to a reconstructed magnetic 
field map in an idealized case with no edge removal in the absence of numerical errors. 

4. Working with the field ;Brec(7), following the steps outlined in Sec. IIV Al and Sec. IIVBI we 
build the unbiased estimator P^^ and calculate its covariance matrix. The resulting estima- 
tor is equivalent to one that could be constructed in an ideal situation, without numerical 
errors, in which we could have worked without edge removal. 
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FIG. 10: The averaged values (red larger dots) and error bars (red larger error bars) following from 
simulations of the unbiased estimators of the power spectrum i(i. + l)C^'^(total)/27r (left panel) and 
iij + l)C^^(gw)/27r (right panel). In both panels, the black solid line denotes the theoretical values of 
the underlying power spectra. For comparison, in both panels, we plot the averaged values (blue smaller 
dots) and simulated error bars (blue smaller error bars) of the unbiased estimators for an ideal case without 
information loss (see text for the details). In both the panels, we have considered a case with no instrumental 
noise. 



The resulting averaged value for the estimator (blue dots) and the corresponding error bars (blue 
error bars) are plotted in Fig. [TOl Once again, we find that the average values of the estimators 
are practically coincident with the theoretical (input) values. In both the panels, the blue error 
bars are slightly smaller than the red ones for all multipole bins. The difference reflects the loss of 
information due to edge removal. In Fig. [TTl we plot the ratio of the two error bars as a function 
of the multipole bin. This ratio in almost everywhere less than 1.2. The signal-to-noise ratio (j48p 
calculated for the ideal is SjN = 9.59, which is less than 15% higher than the practically relevant 
example considered in the previous subsection. The results of this section demonstrate that the 
loss of information, gauged by the increase in the error bars of the spectral estimator, is sufficiently 
smaU < 15%. 
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FIG. 11: The ratio of practically achievable error bars (red error bars in Fig. [TU|) to the corresponding error 
bars in an information lossless case (blue error bars in Fig. llOp , as a function of multipole £. 

D. Power spectrum estimators in the presence of instrumental noise 

In the previous subsections we have considered a situation in which the magnetic type of polar- 
ization was generated solely by gravitational waves and cosmic lensing. In realistic observations, 
in addition to these two contributions, there are contaminating contributions from various other 
sources like instrumental noise and astrophysical foregrounds. This said, it is reasonable to assume 
that, for an appropriate choice of observed sky region, the astrophysical foregrounds are typically 
expected to be small in comparison with instrumental noises [5]. For this reason, we shall ignore 
the foreground contaminations, and restrict our analysis to the study of power spectrum estimators 
in the presence of only instrumental noises. 

The pseudo estimator Di (|37p has the following expectation in the presence of noise (compare 
with no noise case (l38l)') 



(De) = J2 Mu.NJBlC^'' + (A/f ^) , (49) 

where MP is the pseudo estimator for the full sky noise power spectrum Nf. The expectation 
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value of this noise estimator is 

The presence of noise leads to a redefinition of the unbiased estimator P^ 

Pr = E M^JY^P^^eiDe - (A/f ^)), (51) 

b' e 

with matrices pf,ii and M^fe/ given in (j44p and (|45p . The covariance matrix for this estimator has 

the form given by (j47p . where (ADiADci) are calculated from the right side of (j40p with BfCf 

terms replaced by {BJCf^ + N^^y 

The estimator P^ defined in (15ip is an unbiased estimator for the i?-mode of polarization 

power spectrum £(^+ l)Cf /2tt, where Cf contains contribution from both gravitational waves 

(gw) and cosmic lensing (lens) 

C7fS = C7f^(gw) + Cr(lens). 

However, if we are primarily interested in detection of gravitational waves, we can treat the cosmic 
lensing contribution as effective noise, and define an unbiased estimator for the S-mode power 
spectrum due to gravitational waves (.{(. + l)Cf^{gw)/27r as 

Pf ^(gw) = ^ M-,}Y,pMDe - {MD), (52) 

b' i 

where the effective noise term {M^^) contains contribution from instrumental noises and cosmic 
lensing 

(AT/^) = ^ Mu'NJ (S|,C#^(len) + N^'^). (53) 

I' 

The covariance matrix for this estimator is same as that calculated for estimator (|5ip . 

E. Expected performance of ground-based CMB experiments 

In this subsection we shall investigate the prospects of detecting the i?-mode signature from 
relic gravitational waves by two future ground based experiments, QUIET [14] and POLARBEAR 

Let us firstly consider the QUIET experiment. We shall restrict our analysis to the 40GHz 
frequency channel. The FWHM for the Gaussian beam at this channel is Op = 23', and the 
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60(1). We shall assume that 



expected instrumental noise is N^ = 2.72 x 10~^^K^ [ij] (see also 
experiment will observe /sky = 3% fraction of the sky, corresponding to ^o = 20°. following 
the steps outlined in Sec. IIV Al and Sec. IIVB|, using the experimental characteristics for QUIET 
experiment, we construct the unbiased estimators P^^ and P^'^(gw) and their covariance matrices 
for 1000 realizations with r = 0.1 and r = 0.01. The average values for the estimators and their 
corresponding error bars are plotted in Fig. W2\ for r = 0.1 (left panel) and r = 0.01 (right 
panel). The error bars in this case are larger than the error bars in Fig. [10] due to the inclusion of 
instrumental noises. The total signal-to- noise ratio in (j48p is S/N = 7.05 for r = 0.1 model, and 
S/N = 1.25 for the model with r = 0.01. 

We now turn to the POLARBEAR experiment. Once again, we restrict our study to the 
performance of the best frequency channel at 150GHz. The FWHM for the Gaussian beam is 



m- 



As above, we assume 



Op = 4', and the expected instrumental noise is N^^ = 4.22 x 10~^^K^ 
/sky = 3%. In order to study the performance of POLARBEAR, we simulate 1000 realizations of (Q, 
U) maps with r = 0.1. Before proceeding to construct the power spectrum estimators, one should 
notice that, in comparison with QUIET, the value of Op = 4' for POLARBEAR is substantially 
smaller. Thus, in order to avoid leakage from higher multipole electric type polarization, we 
firstly apply the 'map smoothing' procedure outline in Appendix |Dj Following this, we construct 
the estimators P^^ and P^^(gw) and their covariance matrices following the steps explained in 
Sec. IIV Al and Sec. IIVBI In Fig. [13] we plot the average values of the estimators and their error 
bars. The error bars for POLARBEAR experiment are considerably larger than those in Fig. [12] 
(and Fig. [TO]l due to larger instrumental noise in comparison with QUIET. Finally, we calculate 
the signal-to-noise for the POLARBEAR experiment to be S/N = 4.31 for a model with r = 0.1. 
It is worth pointing out that, although in our estimations above we relied on the performances 
of a single best frequency channel for QUIET and POLARBEAR, these experiments will observe 
in several frequency channels. The combining of data from several frequency channels will have 
an effect of reducing the total effective instrumental noise. In addition, these experiments could 
potentially observe larger portions of sky. Both these points could potentially increase the detection 
ability of these experiments. On the other hand, one should remember that various foregrounds 



63| and systematic errors 6j] would increase the effective noise, thereby reducing the detection 
ability. One should remember these caveats, when looking at various signal-to- noise estimates, 
including the ones presented above. 

At the end of this subsection we shall briefly discuss a widely used analytical approximation 



for signal-to-noise. In this approximation S/N ex y^fsky, where /sky is the sky cut factor. This 
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approximation follows from following considerations. In the case of full sky coverage, one can 
construct an unbiased estimator Df-^ (where X = T,E or B) for the various power spectra C^-^ 
in a straightforward manner (see for example 



35|, 



65| for details) . In this case the covariance matrix 



is diagonal with 



{^Df^^Df^) = \l :^^{Cr + N^ Bf) 



with {21 -\- 1) in the denominator on the right side playing the role of number of degrees of free- 
dom for a given multipole L The above expression was elegantly extrapolated for temperature 
anisotropy power spectrum estimator DJ to partial sky surveys in 66|. The author proposed to 



replace (2^ + 1) with the effective number of degrees of freedom (2^ + l)/sky in the above expres- 
sion, to account for the loss of information that arrises due to partial sky coverage. This simple 
consideration was extended to i?-mode power spectrum estimator in [671], and was further extended 
to account for multipole binning [68]. These approximation lead to 



with i being the central multipole in each bin. In this approximation, the total signal-to-noise ratio 
for gravitational wave signal in the i3-mode of polarization takes the form 



S/N 



y^ Mpngw)) y 



In order to guage the performance of this approximation, in Fig. [TT] and Fig. [T21 we plot the 
error bars calculated using (j54p (grey error bars). For this calculation we have set /g^y = 0.024 
corresponding to an effective top-hat window function with 6q = 18°. One can see that, in both the 
figures, the analytical approximation leads to smaller error bars than those obtained from numerical 
simulations. We use (|55p to calculate the analytical signal-to-noise ratio for the two considered 
experiments. The results for signal-to-noise ratio are summarized in Table HI It can be seen that 

y overestimates the detection 



39, 



69[ have pointed out that 



the analytical approximation for signal-to-noise ratio (j55p considerabl 

ability, particularly for smaller values of actual S/N . Several works 

the analytical approximation ()55p exaggerates the detection ability. However, these paper argued 

that the primary reason for overestimation is due to the omission of possible contaminations from 

£^i3-mixing. However, our approach shows that the analytical approximation (I55p with an effective 

sky-cut factor also overvalues the S/N in comparison with the case with no EB-raixhig. One should 

therefore use this approximation with caution [75]. At the same time, it is very important to point 
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FIG. 12: The averaged values and simulated error bars of the unbiased estimators for the power spectrum 
£{£+ l)Cf^ /2tt (green dots and green largers error bars) and £{£+ l)Cf^{gw)/2TT (red dots and red larger 
error bars). In both panels, the solid lines denote the theoretical values for these power spectra. In this 
figure, we have considered the instrumental noise for QUIET experiment. The left panel shows the results 
for an input cosmological model with r — 0.1, while the right panel shows the results for an r = 0.01 model. 
In both panels, the smaller error bars calculated using the analytical approximation (j54p are plotted in grey. 



out that this conclusion about overestimation is based on analysis of small sky coverage and the 
use of pseudo-Cf estimators with the uniform weight function VV(7). In contrast, for large scale 



surveys 



54( 1 or the small scale surveys by using the pseudo-C^ estimators with the optimal choice of 



the weight function W{'j) [41], the conclusion might change. Especially, for the large scale surveys 
and maximum likelihood estimators, the discussed analytical approximation may underestimate 
the true S/N, as was shown for temperature anisotropy in [54| . 
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FIG. 13: The results for POLARBEAR experiment. The averaged values and simulated error bars of the 
unbiased estimators for the power spectrum l{i + l)Cf^ /2tt (green dots and green larger error bars) and 
£{£ + l)Cf^{gw)/2TT (red dots and red larger error bars), calculated for an input model with r = 0.1. The 
solid lines denote the theoretical values for these power spectra. The smaller error bars calculated using the 
analytical approximation ([54)) are plotted in grey. 



TABLE I: The total signal-to-noise S/N for the gravitational waves signal in the B-mode of polarization 
for the various cases considered in the text 





ideal no noise case 
r = 0.1 


QUIET noise 
r = 0.1 


QUIET noise 
r = G.01 


POLARBEAR noise 
r = 0.1 


simulated S/N 


8.26 


7.05 


1.25 


4.31 


analytical S/N 


11.24 


10.76 


3.33 


7.15 



V. CONCLUSION 

In this paper we have proposed a new method to construct pure electric and magnetic type 
fields <?(7) and ^(7) from polarization field given on an incomplete sky. Due to the differential 
definitions of these fields, we avoid the so-called ii^i^-mixing problem. In practice when working 
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with pixelized maps, residual leakages from various numerical errors require the removal of data 
from a narrow edge on the boundary of the observed sky. This leads to a minor loss of information 
in comparison with the idealized lossless case considered in Sec. llVCi 

A major advantage of our approach is that the constructed fields £{"y) and ^(7) are scalar. For 
this reason, the various techniques developed for the analysis of temperature anisotropy maps can 
be directly applied to these fields. As an important and motivating application, we discuss the 
construction of an unbiased estimator for the i?-mode power spectrum Cf^, using the pseudo- 
Ci estimator approach. We find that our method is computationally feasible even in the case of 
high resolution maps. In particular, it takes 2.5 minutes on a laptop (2.4GHz processor and 2GB 
memory) to perform all of the calculations, including the the calculation of ^8(7) in pixel space 
with A'^^ide = 512, and the construction of unbiased estimators for Cf^. 

With the help of the constructed unbiased estimator, we have investigated the ability to detect 
gravitational waves through the S-mode of polarization in CMB experiment covering 3% of the sky. 
In the absence of instrumental noise, we find S/N = 8.26 for a model with r = 0.1. This value is 14% 
smaller than an idealized situation with no information loss. In the case of realistic experiments, 
the signal to noise reduces to S/N = 7.05 for QUIET and S/N = 4.31 for POLARBEAR. 

In conclusion, we would like to point out that, a similar analysis can be applied to large sky 
surveys. In particular, for Planck satellite and the planned CMBPol experiment, one can construct 
unbiased estimators for the polarization power spectra Cf^ and Cf^, by synthesizing the approach 
outlined in this paper together with the hybrid estimator method suggested in [Sjj. We leave this 
task for future work. 
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Appendix A: Numerical calculation of the correction term ct in pixel space 



In order to calculate the correction term ct in (|'24|) and (|25p . one needs to be calculate the terms 
{QW)x, {QW)y, {UW)x and {UW)y. Below, we discuss the calculation of {QW)x- The other 
terms are calculated in a similar manner. 

We expand the polarization fields {Q+iU)W and {Q—iU)W in terms of spin- weighted harmonics 

(Q(7) ± iUmW{^) = Y, a±2/m ±2Yim{l). (Al) 

em 

It follows that 

QW{^) = - ^ Ee^Xi^erriil) + iBhaX2/'m{l), (A2) 



Im 



where 



Elm = —{a.2/m. + a~2,fo«)/2, Bira = —{a2,lm — a_2/m)/2«, 
Xl^lm = ( 2^m + -2^m)/2, X2/m = ( 2^m — -2^m)/2. 

The quantity {QW)x can be numerically calculated using 

{QW)x = diQW)/d9 = -Yl Ee^{dXi^e^/de) + iBem{dX2,em/d9). (A3) 

£m 

Thus, using the expansion coefficients in (jAip and expression ()A3P one can calculate the quantity 
{QW)x in terms of quantities QW, UW and functions [dXn^im/dO). We would like mention here 
that in the HEALPix version 1.23, the subroutine alm-map_template.f90 had a bug, that led to 
erroneous results for {QW)x and {UW)x [72i] . This problem has been fixed in the latest HEALPix 
version 1.24. 

In the present paper we use a simple analytical form (j27p for the mask window function W{'y). 
For this window function, the various derivatives, such as Wxi'j) and Wxxil), can be calculated 
analytically. However, in practical situations, the window function does not have such a simple 
form (see for instant [73]). For this reason, one would need to calculate the various derivative 
terms, Wx, Wy, Wxx, Wyy and Wxy, numerically. This can be done in the following way. One 
firstly defines the multiple expansion coefficients Wim in the standard way 



W,m^Jw{j)YUj)d^. 
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Following this, one calculates 




dW -5-^ 






{ ^ 



The other quantities can be calculated in an anologous manner. It is important to point out that 
the steps mentioned above can be realized in a straightforward manner using the anafast and 
synfast routines in the HEALPix package. 

Appendix B: Construction of magnetic map B^odj) from simulated polarization maps 

In this appendix, we outline the steps which were used to simulate the polarization maps and 
construct the pure magnetic map ;Srec(7)- 

1. We generate the mask window function 1^(7) using (|27p in pixel space using the standard 
pixelization scheme used in HEALPix with Ngi^e = 512 (or Nside = 1024 in the example in 
SeciniB]). 

2. Using synfast HEALPix routine, we generate full sky ((5(7), U{'y)) maps with A^side = 
512 or 1024, using cosmological parameters pO]) and appropriate value of tensor-to-scalar 
ratio r as input. Using the window function ^^(7), we build the masked (Qij), U{'y)) maps 
(where Q = QW and U = UW). 

3. Using anafast HEALPix routine, we calculate the coefficients (-E^m; -B^m)- The field B{^) is 
calculated from B^j^ according to p^ using synfast routine. 



4. With coefficients (-E^m, -B^m) we construct the fields QW{^), UW{^), {UW)x{'^) and 
{QW)y{'^), using the 5*^^ option in synfast routine. 

5. Using the fields UW{^), {UW)x{'^) and {QW)y{'^) constructed in the previous step and 
analytical expressions for W, Wx and Wxx (derived by differentiating ([27|) ). we calculate 
Im(ct) in ([29)1 . 



6. The pure magnetic field Brecil) is now constructed from B, W and Im(ct) using ([32]) . The 
pure magnetic field Bredl) is truncated at the edges in order to remove residual leakages 
associated with numerical errors. 
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Appendix C: Pseudo Estimators for a special choice of weight function yV(7) = W^(7) 

In Sec. IIV Al it was pointed out that, in principle one can construct pseudo estimators of the 
power spectrum by adopting an arbitrary weight function W{'j) in (|33p . Above, in the main text, 
we have focused on a specific case corresponding to a uniform weight function W(7) = 1 . This 



choice is near 
multipoles 



y optimal for high multipoles. However, this choice becomes sub-optimal at lower 
4ll |. In this appendix we study another possible choice for the weight function, namely 
W{^) = VF(7), where W{^) is the mask window function in ([27D with ^o = 20° and 6i = 10°. Note 
that the function 1^(7) is the same window function that was used for constructing Brcci'l)- With 
this choice, the resulting pseudo estimator is equivalent to the pure B-mode estimator studied in 

Q. 

The construction of the pseudo estimators and the corresponding unbiased estimators follows 
closely the discussion in Sec. llVi The only difference is that the definition of coefficients aim in 
(l34ll are modified to 



aim = j d^ BrUlWiWLil)^ (CI) 

and the quantities w'{^) and w'^ in ([36]) and ([39|) would now be replaced by W and its power 
spectrum, respectively. 

In Fig. [TH we plot the unbiased estimators for the power spectra i{C + l)C^^(total)/27r (left 
panel) and £(^-|-l)C^^(gw)/27r (right panel) together with the corresponding error bars (thin blue 
error bars). It can be seen that, in comparison with the estimators in the case of a uniform weight 
function, the error bars of the new estimators are larger at high multipoles, but are smaller at low 

4l[ | that the optimal weight functions for high 



multipoles. This result is consistent with findings in 

multipoles tend to the top-hat function. On the other hand, for low multipoles, the optimal weight 

function tend to smooth out (see the right panel of Fig. 2 in 4l| for a concrete example). 



Appendix D: Smoothing the polarization maps 

The high value of power spectrum Cf^ of the electric component at large values of multipoles 
(due to the presence of Nf factor) leads to substantial leakage of power into the reconstructed pure 
magnetic field ^Srec- This leakage seeps through to low multipoles playing a role of residual effective 
noise. In order to reduce this contamination, below we introduce a map smoothing procedure for 
polarization maps. The idea behind this method is similar to the 'prewhitening' method suggested 
in y. 
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FIG. 14: The averaged values (blue smaller and red larger dots) and error bars (blue smaller and red larger er- 
ror bars) following from simulations of the unbiased estimators of the power spectrum £{£ + l)Cf ^ (total) /2tt 
(left panel) and £{i+l)Cf^{gw)/2Tr (right panel). In both panels, the black solid line denotes the theoretical 
values of the underlying power spectra. The blue smaller dots and error bars denote the result by adopting 
the weight function W(7) = W{'j), and the red larger dots and error bars denote the result by adopting a 
uniform weight function W(7). Note that, in both panels, we have considered a case with no instrumental 
noise. The red larger dots and error bars are identical to those in Fig. [TU] 

The smoothing procedure is simple and straightforward in the case of full sky coverage. Given 
the (Q, U) polarization maps on a full sky, one can calculate the multipole coefficients Eirn and 
Bum using (dJ-dH). In order to smooth the polarization maps we firstly use a damping function to 
modify the multipole coefficients 



E' 



Im 






B' 



£m 






(Dl) 



Following this, we reconstruct the smoothed polarization fields Q' and U' using the modified 
coefficients E'^^ and B'^^ in the standard way. The reconstructed maps can be thought of as the 
the result of observing the original (Q, U) polarization field in an experiment with FWHM of the 
Gaussian beam equal to Op- Overall, the smoothing has effect of exponentially damping the power 
in high multipoles. 
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We can extend the smoothing procedure to the case of partial sky coverage. Given the {Q, U) 
polarization maps on a partial sky, we calculate the coefficients E^j^ and i?£m- These coefficients 
are smoothed in analogy with full sky case 



i^^"! ~ ~ _it^el 



E'lrn = Elme 2 8 in 2, 5'£„ = 5fo„e 2 8 In 2. (D2) 

The smoothed polarization maps (Q', U') are reconstructed from the modified multipole coefficients 
E'irn and B' i^a- 

It is important to point out that, in the case of partial sky coverage, the smoothing procedure 
outlined above introduces mixture of electric and magnetic polarizations. In particular, even if the 
original {Q, U) did not contain magnetic type of polarization, the smoothed map {Q' , U') would 
contain it. We have verified numerically that in practically interesting cases the resulting mixture 
is very small, and would not significantly affect the ability to detect gravitational waves. 

In order to verify the small of the resulting mixing, we have performed the following calculation. 
Using an input model with no i?-mode of polarization (i.e. Cf^ = 0) and Op = 30' we generated a 
full sky (Q, U) map. Following this, we truncate the map to keep the data from only the northern 
hemisphere. Using the procedure outlined in Sec. lUand Sec. IIIII we construct the field B^ec-, which 
is expected to be equal to zero except for the residual leakage. The psedo power spectrum D^ for 
this field (red line) is plotted in Fig. [151 Following this, for the same input model, we generate the 
(Q, U) map with Op = 10', once again restricting the data to just the northern hemisphere. We 
now smooth this map with Op = 30' using the anafast, alteralm and synfast HEALPix routines. We 
construct B'^.^^ from the smoothed {Q' , U') map and plot the corresponding pseudo power spectrum 
D'^ (green line) in Fig. [T5j The difference between the two spectra D£ and D'^ can be interpreted 
as the result of mixing introduced by map smoothing (blue line in Fig. [T5|) . It can be seen that, 
the mixing due to smoothing is QUIET small, comparable to leakage due to numerical errors, at 
all the relevant multipoles. The power spectrum of the leakage due to smoothing peaks at £ ~ 400. 
It can be completely neglected at low multipoles {£ < 50). 

Below we give a heuristic argument to understand these results. The two sets of multipole 
coefficients {Eim, B^^) and {Eim, B^^) are related by the following expression (see for instant 

0) 

Elm + iB^m = 2_^ 2l(im)(i'm')[E&m' + iBi/m'], (D3) 



i'm' 



where 2l(£m){e'm') is the coupling matrix, which depends only on the mask window function. From 
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this relation, it formally follows that 



i'm' 

For the smoothed multipole coefficients one has 



E'^^ + iB'g^ = 2^(2-^ ){em){e'm')[E'e'm' +iB'i' 



I'm' 



ll^2lk 



= 2^(2-^ ^){e.m){em')[Eem' +'i'Bii^rn']e 2 8ln2 . 
i'm' 

,n 

Since the coupling matrix 2l{em){e'm') is sharply peaked at £ = r [32], the above expression can be 
approximated by 

^-^ ~ ~ 1 *^*F 

B'em + 'i'Bem ^ 2^(2^" )(^m)(rm')[-^^'m' + «-B£',m']e~2 8hr2 

i'm' 

It therefore follows that 

E'lm ~ Eime~^-sT^, B'f,^ ^ Berne'^^^. (D5) 



I2„2_ 



1 i'^e-p 
Thus, in this approximation, there is no mixing introduced by smoothing. Moreover, e 2 8 in 2 ss 1 

for £' <C l/Op, leading to E'^^^ Ri Eim, B'^^ ss Bim for low multipoles. In reality the coupling matrix 

2l(im)(i'm') is often non-invertible, and therefore ()D4p might not be valid. However, the resulting 

leakage is still small, and completely negligible for I <C ^/Of- 
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(i.e. line 1) shows the pseudo power spectrum Df^. Note that, the black line (i.e. line 1) and red line (i.e. 
line 2) are identical to those in Fig. |3l 
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